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Abstract 

A new statistical model is proposed for the geomagnetic secular variation over the past 5Ma. 
Unlike previous models, which have been based on the assumption of a Fisher distribution for 
either VGP positions and/or field directions, the model makes use of statistical characteristics of 
the present day geomagnetic field. The spatial power spectrum of the non-dipole field is con- 
sistent with a white source near the core-mantle boundary with Gaussian distribution. After a 
suitable scaling, the spherical harmonic coefiicients may be regarded as statistical samples from a 
single giant Gaussian process; this is our model of the non-dipole field. We can combine this 
model with an arbitrary statistical description of the dipole and compute the probability density 
functions and cumulative distribution functions for declination and inclination that would be 
observed at any site on the surface of the Earth. Global paleomagnetic data spanning the past 
5Ma are used to constrain the statistics of the dipole part of the field. A simple model is found to 
be consistent with the available data: (a) with two exceptions, each Gauss coefficient is indepen- 
dently normally distributed with zero mean and standard deviation for the non-dipole terms com- 
mensurate with a white source at the core surface; (b) the exceptions are the axial dipole, g^ , and 
axial quadrupole, , terms; the axial dipole distribution is bimodal and symmetric, resembling a 
combination of two normal distributions with centers close to the present-day value, and its 

sign-reversed counterpart; (c) the standard deviations of the non-axial dipole terms, gl and A/, 
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and of the magnitude of the axial dipole are all about 10% of the present-day component; (d) 
the axial quadrupole reverses sign with the axial dipole and has a mean magnitude of 6% of its 
mean. An advantage of specifying the model in terms of the spherical harmonic coefficients is 
that it is a complete statistical description of the geomagnetic field, enabling us to test specific 
properties for a general description; both intensity and directional data distributions may be 
tested to see if they satisfy the expected model distributions. 

Introduction 

It is by now generally accepted that the geomagnetic field is generated by some sort of 
dynamo process in the Earth’s core (see e,g., Moffatt, 1978). The field is by no means constant; 
changes on time scales from ten years to ten thousand years are classified as secular variations. 
Such fluctuations are consequences of the nonsteady nature of the field-producing mechanism 
within the core. More rapid variations, if they exist, cannot be detected because of screening due 
to mantle conductivity and interference from fluctuations of external origin. At the longer 
periods, it is difficult to dissociate mere secular variation from the actual reversal of the main 
field and it is quite possible the distinction is artificial. Unfortunately, the theory of the geo- 
dynamo is not currently in an advanced enough state to make more than the simplest predictions 
about the paleofield variation {e.g., Merrill and McElhinny, 1983, Chapter 9). When a system is 
as complex in space and time as the geomagnetic field, it is a natural and efficient strategy to call 
upon some kind of statistical characterization of it to uncover general properties, rather than to 
attempt a description of the behavior of the system at every instant and at all points. Previous 
statistical descriptions of the paleomagnetic field have centered around the use of Fisher statistics 
(Fisher, 1953) and attempts to describe the observed variation in geomagnetic field dispersion 
(from paleomagnetic data) as a function of latitude (Creer et aL 1959; Creer 1962; Cox 1962, 
1970; Baag & Helsley 1974; McElhinny & Merrill 1975; Harrison 1980; McFadden & McElhinny 
1984). The most successful of these to date has been the so-called Model F of McFadden & 
McElhinny (1984), for which the average intensity of non-dipole components of the field at any 
latitude are assumed proportional to that of the dipole field, and the angular distribution of VGP 


-3- 




directions is assumed to exhibit axial symmetry. The model provides a good fit to the variation 
in VGP angular dispersion as a function of latitude. The aim of tUs work is to provide a statisti- 
cal description of the field variation, by combining the properties of the present-day field and of 
paleomagnetic measurements. This statistical description does not just describe the angular 
dispersibn expected at any latitude, but provides complete statistical distributions expected for 
the secular variation at any site on the Earth’s surface. The resulting models are specified in 
terms of the spherical harmonic coefficients describing the geomagnetic field. 

The general plan of this paper is as follows. Two kinds of data describing the main geomag- 
netic field each supply components of the model and these are discussed; they are the modern 
data, with their excellent geographical coverage of the Earth, and the paleosecular variation data 
which span long time periods, but in general have poor global coverage. The modern data pro- 
vide a good estimate of the spatial power spectrum of the geomagnetic field, stnd enable us to 
show that the non-dipole part of the field is consistent with a white source near the core-mantle 
boundary. By suitable scaling of the spherical harmonic coefiicients it is possible to regard them 
as statistical samples from a single giant Gaussian distribution. Supposing that a kind of unifor- 
mitarianism has held true, we assume that this giant Gaussian process has always, been a good 
description for the non-dipole geomagnetic field. However, the modern data cannot help us with 
the statistical variation of the dipole, which involves much longer time spans; we must use psJeo- 
data to study its behavior. Using a statistical framework it is possible to predict the cumulative 
distribution functions for quantities related to inclination and declination, which are the most 
accurately measured paleodata. The paleodata provide us with empirical cumulative distribution 
functions for the directional data from Hawaii and Iceland, where there are extensive paleosecular 
variation data from lava flows, and also for a global data set based on a compilation by Lee 
(1983). We construct a simple model with a number of free parameters, and adjust them to 
reconcile the model and the paleodata. 
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1. Short-Term Data and Nondipole Field Statistics 

Satellite and observatory measurements and historical records of the geomagnetic field fall 
into the category of short-term secular variation data. Observatory data rarely extend back more 
than 100 years, but provide reasonable simultaneous global coverage for that time (although the 
distribution of sites is heavily biased towards the northern hemisphere and Europe, and per- 
manent observatories are necessarily restricted to land which is very unevenly distributed over 
the surface of the Earth ). Satellite measurements have been performed at irregular intervals 
since the launch of Sputnik 3 in 1958 (see Langel, 1985 for a complete list), and the resulting data 
have varied in quality and extent of global coverage. Magsat (Langel et ai, 1980, 1982) provided 
the first truly global survey of the vector components of the geomagnetic field. 

The method of analysis for these data was devised in 1838 by Gauss and involves represent- 
ing the geomagnetic field of internal origin as the gradient of a potential 4^ whose spatial varia- 
tion 'is specified by an infinite sum of spherical harmonics. 

and are known as the Schmidt partially normalized Gauss coefficients and provide a com- 
plete description of the geomagnetic field. Appendix 1 gives the relationships between fully and 
partially normalized spherical harmonic representations and indicates how they are computed. In 
this work we make use of the GSFC980 spherical harmonic model of Langel et al. (1980), derived 
from Magsat data. 

Previous secular variation studies have been seriously hampered by an inability to handle 
the non-dipole components of the field in a quantitative manner. We propose a way out of the 
difficulty based upon a simple observation of the present-day field, that the non-dipole terms can 
be described by a zero-mean Gaussian process, and that it is plausible to assume this property 
has held for all time. 

Mauersberger (1956) first defined a kind of geomagnetic power spectrum based on the spheri- 
cal harmonic expansion of the geomagnetic field. Subsequently, other workers have used various 
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forms for the spectrum with different weighting as a function of degree. Here the form devised by 
Lowes (1974) is used. For each spherical harmonic of degree /, the power at radius r is Ri, 

Ri (r)= (^)»C^»)(/ + 1) 2 { (gr? + [hr? } 

m-O 

where gr hr ^be Gauss coefficients in the Schmidt quasi-normalized spherical harmonic 
expansion. Equivalently, 

Ri = <B i‘B I > 

where Bt is the magnetic field associated with degree I in the expansion and < > denotes the 
average over the surface of the Earth. 

At the Earth’s surface the spectrum drops exponentially with a slope indicating a white 
source approximately at the core’s surface (see Figure 1). Lowes gives: 

R, ~ 4.0X10® (4.5)-'(|xr)» 

The deviations from this law for I > 8 are thought to represent power from crustal magnetiza- 
tion. The spectrum has subsequently been refined, through the use of the more complete Magsat 
data set (Langel & Estes 1982). They estimate that the spectrum is consistent with a white source 
174 km below the seismic core-mantle boundary. 

If the spectrum of the geomagnetic field were precisely the same throughout geologic time 
then secular variation would just be reflected by a rearrangement within the individual gr and 
hr while the power Ri remained constant. This implies a kind of 2/ + 1-dimensional Fisherian dis- 
tribution for the Gauss coefificients. However, exactly constant power in each spherical harmonic 
degree contradicts observation; continual small changes in Rt are undoubtedly occurring. We 
therefore choose a description in which the power is also a random variable. The simplest model 
treats the individual spherical harmonic coefficients as normally distributed random variables. 
This simple model or a similar assumption has been used by other workers as a statistical descrip- 
tion of the secular variation (see e.g., Gubbins, 1983, Eckhardt, 1984) and has some very attrac- 
tive features. For example, the statistical prescription is invariant under reorientation of the 
coordinate axes. Many other specifications force a privileged status on one particular axis system. 
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something that is clearly undesirable. Another attractive feature is that the model is amenable to 
testing, at least for the present day field. 

We use the Lowes spectrum as a guide in construction of our statistical model. For the 
non-dipole terms we choose a model spectrum that is exactly flat at the surface of the core: 

E{i2J = E{ 2 (/+1) + (*/”)“)> 

m “0 

= / S: 2 

where E{ } is the expectation of the parameter in the brackets, c/a is the ratio of the core radius 
to that of the Earth (0.547) and a=27.7|xT is a fitted parameter. In addition we assume that 
within each degree /, gf* and are independent identically distributed normal random variables 
of zero mean. If var(A|”*) = var(( 7 |”*) = <r * then 

2 _ (c/g 

“ (/ + 1)(2/ + 1) 

If we define scaled Gauss coefficients 

§r = yi^r 
hr = Yihr 

where V|^=(/ + l)(2/ + l)(a/c)*^, then these coefficients (which completely characterize the non- 
dipole field for statistical purposes) are just independent samples of a single zero-mean Gaussian 
process with variance a*. Even though there is some evidence to indicate that the level at which 
the spectrum is white lies somewhat below the core mantle boundary we prefer to characterize 
the spectrum in terms of the known physical constants c and a . The reason for this is that we do 
not believe that the true level is very well constrained. In any case, the method used here in 
fitting the data yields an excellent fit to both the spectrum (see the straight line fit on Figure 1) 
and, as we will see, the Gaussian distributional form for the non-dipole field. Subsequently, we 
may want to allow the level where the spectrum appears white to vary. This may be done by 
changing the variance of the scaled Gaussian distribution. 
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The spherical harmonic field coefficients for 1980 (GSFC980) data for degree / =* 2 through 
to / = 8 or higher may be used to test whether this Gaussian distribution is actually observed in 
the appropriately scaled coefficients. Figure 2 shows the cumulative distribution versus the 
theoretical normal curve. The agreement between the two curves is impressive. While there is no 
m)§ans‘of determining a unique statistical model from so few data, it is possible to test whether a 
data sample is consistent with an assumed distribution. The Kolmogorov-Smirnov test (see e.g. 
Massey 1951, Kendall & Stuart, 1979) is a powerful means of determining probabilistic bounds 
within which samples from a particular theoretical distribution function should lie; one virtue of 
it is that the test itself is independent of the underlying theoretical distribution function 
(although of course this must be known to compute the discrepancy statistic djy). The data of 
Figure 2 are consistent with an underlying normal distribution for the non-dipole field; they exhi- 
bit a maximum discrepancy from the theoretical normal distribution of i^^=.061, where 7V=77- 
However, if the three dipole terms associated with / = 1 are included, we find the K-S test for 
normality fails, demonstrating that the dipole terms do not follow the same pattern as the rest 

of the field. One might suspect that just the axial {ue., gi) part of the dipole field is anomalous 
and we can test this by including the non- axial dipole terms and hi) in the tested distribu- 
tion. Then the population tests out as Gaussian at the same level as the non-dipole field. There- 
fore it appears that according to this model, only the axial dipole term gi is anomalous. 

As a description of the non-dipole field, our model has simplicity and economy; only one 
number, the variance, needs to be determined empirically. It seems extremely unlikely that such 
a simple description is an accident and pertains to today’s field alone. Rather, it is a reasonable 
working hypothesis that this model accurately provides the statistics of the non-dipole part of the 
paleofield as well. Now, for example, we can compute the statistics of the ”noise” contributed to 
paleomagnetic measurements due to non-dipole sources. F rom the spherical harmonic representa- 
tion we can show 1) that the components of the non-dipole field at the Earth’s surface are also 
independent Gaussian variables with zero means, 2) that the variances of the horizontal com- 
ponents locally are equal, and 3) that they are different from that of the vertical component. 
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This is done in Appendix 2. The inequality of the vertical and horizontal variances makes it evi* 
dent that the distribution of directions for the non-dipole part of the field is not Fisherian. This is 
entirely consistent with observation. 

2. Paleodata and the Dipole Field 

With paleomagnetic data it is never possible to obtain the detailed instantaneous descrip- 
tion of the magnetic field that we have for the present day. However, we can obtain information 
about much longer term variations, which is essential for a description of the dipole part of the 
field variation. The bulk of the data available consists of directional measurements of the 
paleomagnetic field. Measurements of the absolute paleointensity of the Earth’s magnetic field 
can be obtained from lava flows and baked archeomagnetic materials. Such measurements are 
exceedingly time consuming and therefore are far less numerous than directional measurements. 
Nevertheless, they are important if we are to be able to determine average values of field inten- 
sity as well as direction. The available data have been compiled by McElhinny & Senanayake 
(1982) for the past 50 thousand years and by McFadden & McElhinny (1982) for the past 5Ma. 
For directional data, there are two major sources of paleodata that span long time periods, lava 
flows and sedimentary sequences. 

Samples from lava flows typically produce very high quality data, but have the disadvantage 
that they only provide instantaneous recordings of the field properties at the time of acquisition 
of the remanence. This is not, however, an insurmountable problem for a statistical model of the 
geomagnetic field, since we are seeking a probability distribution rather than a deterministic 
model of magnetic field changes with time. Precise age control is not necessary, nor indeed is a 
continuous record of the field; all that is required is that we obtain a representative sample of the 
complete range of the secular variation at a given site. There are two places which could conceiv- 
ably offer sufficient data, namely Hawaii and Iceland. Both islands have experienced repeated vol- 
canic eruptions at relatively short time intervals over the past few million years. Extensive data 
sets are available from these sites already (for Hawaii see Doell and Cox, 1965; Doell, 1969; Doell, 
1972a,b,c; Doell and Dalrymple, 1973; for Iceland see Watkins, McDougall and Kristjansson, 
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1977; Kristjansson tt aL, 1980). From the pool of declination and inclination data for each of 
these sites, we can obtain empirical cumulative distribution functions (cdfs) for declination and 
inclination at the latitude for each of these sites. Figure 3 shows these empirical functions, for 
Hawaii using the data of Doell & Cox (1965), Doell (1972a, b,c) and Doell & Dalrymple (1973), 
and for Iceland using the data of Watkins, McDougall & Kristjansson (1977) and Kristjansson et 
al. (1980). Data from lavas more than 5Ma old were excluded as were those sites with Fisherian 
ags>20°, in order to reduce contamination of the secular variation signal by tectonic movements 
and rock magnetic noise. The inclination data are all mapped onto the lower hemisphere in order 
to avoid the necessity of deciding whether they correspond to a normal or reversed field and the 
declinations are mapped into the range -90®<<l£90°. This could be important for low latitude 
sites, where secular variation might result in an apparent field reversal. The changes of variable 
from declination D and inclination I to their modified forms d and t are performed as follows. 

t = I/I 

and 

d=D for -90“<2>s90* 
d=D + m for -180“<P:S-90° 
d=D-lB0 for 90“</)sl80‘’ 

The other data we have employed comes from the global compilation by Lee (1983) for the 
past 5Ma, which contains about 1100 individual site measurements (as well as many more aver- 
age results not used here), from a mixture of both sedimentary and igneous rocks. In the original 
selection of data for this compilation, data with og6>20'’ were excluded as were those with VGP 
latitudes of less than 45®. This is because low VGP latitudes were considered more likely to 
correspond to anomalous field behavior (Lee, 1983) and the object of the study was to look at 
"normal" secular variation. There are comparatively few sites where there are as many data as 
for Hawaii and Iceland, so generally we have to combine data from a number of different sites and 
obtain cdfs that are averages over a finite latitude band. This will also have the effect of averag- 
ing out non-zonal effects in the secular variation, which might otherwise be a problem even for 
very long records at a single site. Figure 4 shows the empirical cdfs for t, the absolute value of 
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the inclination, when the data are averaged in 8 latitude bands ranging in width from 2° to O**. 
The northern and southern hemispheres are assumed to be symmetric in their behavior so that 
data from equivalent latitudes in the two hemispheres have been combined to yield a better esti- 
mate for the cdfs. This assumption of symmetry, means that any non-zero mean in the odd order 
terms m the spherical harmonic expansion will not be detected. The data distribution is such 
that very few of these bands contain sufficient data to provide an acceptable looking cdf. We 
would expect the cdf to look smooth, and many of these are jagged in shape, especially near the 
tails of the distribution. Two of the bands (11—13° and 52—54°) each only contain data from one 
or two locations and have only 51 and 62 field measurements respectively, and so can hardly be 
expected to sample the whole range of secular variation. In order to quantify the variation in the 
cdfs as a function of latitude we can look at the variation of the expectation and standard devia- 
tion of these distributions. These are computed by numerical integration of the empirical cdfs 
using a Qubic spline quadrature scheme. In the upper part of Figure 5 are plotted the standard 
deviations, o-(d) and o-(t), respectively for d (triangles) and » (squares) for each of the latitude 
bands, whose cdfs are given in Figure 4. In the lower part is plotted the bias, B{t}, in the absolute 
value of the inclination. This is defined as the expected value of the absolute value of the inclina- 
tion as computed from the empirical cdf minus the inclination that would be observed at that 
latitude if the field were due to an axial geocentric dipole. The axial dipole inclination is readily 
computed by 

t'a, = tan~^(2tan\) 

where X is the latitude at the location in question. (Here the average latitude for the band, 
weighted by the number of data at each point, has been used.) Thus the bias is 

B{.} = 

Some idea of the reliability of these distribution parameter estimates may be obtained by using a 
bootstrap technique (of the type reviewed by Efron & Tibshirani, 1986) to estimate one standard 
error in the parameters. The bootstrap method supposes that the distributions of Figure 4 
represent the true underlying distribution of the secular variation data, and finds the standard 
error in the parameter estimate by repeated random resampling from this distribution, computing 
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the parameter estimate for each sample and then finding the sample standard deviation. The 
error bars for the parameters in Figure 5 are one standard error based on 500 bootstrap samples 
from each of the distributions shown in Figure 4. How good these bootstrap estimates of stan- 
dard error are will of course depend on how well the distribution functions of Figure 4 represent 
the*range of secular variation. 

The general picture that emerges for the parameters of the cdfs of Figure 4 is that the incli- 
nation standard deviation rises from about 10° at the equator to around 13° at around a latitude 
of 25° and then slowly decreases with increasing latitude. It should be noted that mapping incli- 
nation onto the lower hemisphere results in a modified inclination standard deviation at the equa- 
tor that is expected to be one half that of the true inclination. The declination standard devia- 
tion increases fairly steadily with latitude as does its uncertainty computed by the bootstrap. 
The solid lines connecting the data on this figure are cubic spline interpolation between data 
points and serve merely to guide the eye. The bists in the inclination is high near the equator 
(this is a consequence of mapping all the inclinations onto the lower hemisphere), but becomes 
consistently negative for latitudes higher than 7 or 8°. The size of the bootstrap error bars sug- 
gests that is would be optimistic to try and interpret the detailed variations in bias with increas- 
ing latitude. 

Both the cdfs with the smallest number of points had somewhat smaller standard deviations 
in comparison with their nearest neighbors in latitude (marked by different symbols in Figure 5), 
and this aroused some suspicion that perhaps the samples were not large enough to contain the 
whole range of secular variation. To check this the data were combined into larger latitude 
bands. In Figure 6 cdfs for the absolute values of inclinations in latitude bands of 15° (for lati- 
tudes 0—45° and a band 25° wide for the highest latitudes) are plotted. These are much smoother 
than the cdfs obtained in Figure 4, and we can feel more secure that they represent an adequate 
sample for the variation. The price paid is that of resolution, the modified inclination standard 
deviations, a(t), computed for these distributions will be increased relative to those obtained 
before, because the expected value for the inclination distribution increases with latitude. The 
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standard deviations and bias for this grouping of the data are shown in Figure 7. They are simi- 
lar in form to those obtmned in Figure 5, but without some of the high frequency variation seen 
there. This supports the idea that some of the features in Figure 5 are an artifact of the small 
number of points used to compute the cdf in some latitude bands. 

The next question that arises is whether this global data set is consistent with those from 
Hawaii and Icelwd. Comparison of the Hawaiian data with that from the corresponding latitude 
band shows them to be quite similar. The global inclination data has a bias of -3.2** compared 
with -4.5® for Hawaii, and the standard deviations are quite similar, 12.0® and 13.2® respectively. 
The same is not true for the Iceland cdf which has a bias of -8.8° and a standard deviation of 
11.8®, compared with a bias of -5.4® and standard deviation of 7.4° for the global data. A look at 
the source of the Icelandic data reveals that there are far more low inclination values than would 
have been expected from the global data set. The cdf for the Icelandic data has risen to a value 
of 0.2 already for inclinations of 60®, while for the global data set this does not occur until past 
85°. This may be a consequence of the fact that in the global data set sites at which the VGP 
latitude lies below 45® are excluded, while we have made no such requirement in our compilation. 
However, inspection of the original data reveals that for the Icelandic data about 10% of the data 
corresponds to VGP latitudes that are less than 45®. If we assume that the lava flows are uni- 
formly distributed in time, and that VGP latitudes this low correspond to times when the field 
was in transition between polarities, then this implies that the field spends 10% of its time in 
transition, or 0.32Ma out of the 3.2 Ma spanned by these two records. There are 15 reversals 
recorded in these lavas, thus if the extrusion rate was uniform, then on average the directional 
part of the magnetic field took about 20 thousand years to complete a reversal. This is 
significantly longer than the average time reversals are believed to take, t.c., between 1 and 
lOthousand years (Merrill & McElhinny 1983, pl64), although Clement & Kent, 1984, suggest 
that there is some evidence that reversals appear to take longer (by a factor of about two) at 
higher latitudes than near the equator. However, the Icelandic cdf could be heavily biaised by 
non-uniform sampling of the field in time due to the erratic nature of the eruptions. 
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A complete description of the statistics of the magnetic field must include time scales long 
enough for reversals, but we hope they will not be so long that the underlying processes must be 
regarded as non-stationary. We choose the past 5Ma as a representative sample for the secular 
variation. This spans about 20 polarity changes (see e.g., Merrill & McElhinny 1983, p.l40), but 
is -not so long that we have to take into account the effects of continental drift. The primary 
paleodata set that we hope to satisfy with our statistical model of the field consists of those data 
from the global data set compiled by Lee (1983) for which individual site measurements are avail- 
able. 


2.1. A Statistical Distribution for the Dipole Field 

From now on we shall assume that a giant Gaussian process of the type described in Section 
1 is an adequate description for the statistical variations in the non-dipole part of the field. The 
question now arises as to what is an appropriate model for the dipole field. This is, in fact, the 
major question addressed in this work. It is not a trivial matter, and indeed, a large part of the 
science of paleomagnetism is devoted to exactly this question: how does the dipole of the 
geomagnetic field behave? As we have been doing all along, we will ignore questions of develop- 
ment with time and concentrate instead on the statistical distribution of the dipole field amongst 
its various possible states. 

Several approaches may be taken to modeling the dipole field variations. These will be dis- 
cussed in order of increasing complexity. In all of the following models (with one exception) it 
will be assumed that the dipole part of the field is statistically independent of the non-dipole 
part. This is clearly not true in general; the complex evolution of secular variation must be con- 
trolled by the geodynamo, in which the dipole and non-dipole variations are presumably inextri- 
cably linked. However, in developing a statistical model for the secular variation we will make 
use of paleomagnetic data, which provides a sporadic sampling of the process throughout time. 
The assumption of independence is then like proposing there exists no requirement for the qua- 
drupole term to always be small when the dipole term is large, or similar relationships between 
other spherical harmonic coefficients. Assuming that the data are not serially correlated 



- 14 - 


(certainly a safer assumption for lava flows than for sediments), the independence assumption 
should not be a serious problem. 

The simplest idea would be to treat the dipole as just another part of the giant Gaussian 
distribution. This is completely incompatible with what is known of field component distribu- 
tions from paleomagnetic data (Cox, 1970 ) in which the axial (t.e., spin-axis-aligned) dipole has 
been dominant throughout geologic time (except during the brief periods of actual reversal). We 
therefore rule out this model. 

The simplest modification of an entirely isotropic dipole model is one that ascribes to the 
dipole variations, the same statistical behavior as the non-dipole field, but, during times not asso- 
ciated with reversal, superimposes a steady axial part to account for the ascendancy of the 
term. This seems a reasonable approach given that the g i and h ^ terms also appeared to satisfy 
the giant Gaussian model for the non-dipole part of the field. Thus the statistical distribution for 
the gi term is envisioned as bimodal and symmetric, a combination of two normal distributions 
with centers at the present day value and its sign reversed counterpart. We take this opportunity 
to introduce some new terminology specifying the magnitude of the zonal terms of the geomag- 
netic field. Let 7° = lj|®l for 1 = 1,2, ■ ■ • ; then the statistical distribution for 7x° will be closely 

approximated by a Gaussian distribution centered on l?°|. (The standard deviation computed 
from the present-day spectrum for the dipole parts of the field is about 20% of the present axial 
dipole magnitude; thus the area under the normal pdf that is truncated by taking the absolute 
value will be very small.) Then at times not associated with reversal the magnetic field com- 
ponents may be written as 

Br' = B, + / 

Be' = 5e + « 

B*' = B* 

where / = 27i°sin\, e = 7 i°cos\ are the mean contributions at a site at latitude X and, as we 
mentioned earlier, B,, B9 and 5*, the random components, are Gaussian, zero mean and indepen- 
dent with variances derived from the spectrum at the core. We can compute the probability 


-15- 


density functions (pdfs) for the commonly measured elements of the geomagnetic field by per- 
forming the necessary integrals. The techniques used for performing these calculations are laid 
out in Appendix 2 for a general form for the axial dipole distribution. Numerical integration of 
those pdfs yield cumulative distribution functions for comparison with the empirical cdfs 
obtained for the data. 

Figure 8 shows the resulting pdfs at a variety of latitudes for t and d, the modified inclina- 
tion and declination. In these calculations we have simply set the mean of the magnitude of the 

axial dipole to today’s value (7° — 30p.T). The distribution for t is skewed at mid and high lati- 
tudes and the standard deviation is large at low latitudes. Figure 9 (solid line) is a plot of the 
standard deviations and the bias as a function of latitude. (The bias is again the expected value 
of 1 minus the axial dipole inclination for that latitude.) Note that even though everything 
except the axial dipole part of the field had zero mean Gauss coefiicients, this does not mean that 
t or even /, the inclination, will average to the axial dipole value. This is fortunate, since there 
exists a considerable body of paleomagnetic evidence in favor of biased inclinations (see e.g. Wil- 
son, 1971; Merrill & McElhinny, 1983 Chapter 6). Also plotted on Figure 9, as symbols, are the 
standard deviations and bias obtained for the global data compilation of Lee (1983), which were 
plotted in Figure 5. It is evident that this model leaves something to be desired. The standard 
deviations are all too high, although the variation with latitude has approximately the right 
shape. The bias in the modified inclination changes sign at too high a latitude (again partly 
reflecting the high standard deviation), and does not reach quite such low values in mid-latitudes 
as it does for the data. 

How closely should we expect the data to follow the model? The error bars in Figure 9 are 
based on the assumption that the empirical distribution functions of Figure 4 represent the true 
distribution functions for the secular variation. The effects of the inadequacy of this assumption 
are hard to measure. Both the t and d data should have slightly higher standard deviations, o-(t) 
and <y(d), than predicted by the model because of the presence of errors in recording and measur- 
ing the field. For the data, a(t) will also be biased towards high values because inclinations from 
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a range of latitudes with slightly differing means are combined to provide a single empirical dis- 
tribution for that latitude band. The only mechanism by which the data standard deviation 
could be lower than that for the model is if there are insufficient data to sample the complete 
range of directions possible for the secular variation. As noted in the section on paleodata, the 
data of Figure 7 suggests that this is only influencing one or two of the standard deviation points 
of Figure 5. These are again marked by asterisks in Figure 9. The net result of the above possible 
sources of bias is that, on the whole we should expect a(t) and a(d) for the model to lie slightly 
below those obtained for the data, except possibly for the data represented by the asterisks. 

From Figure 9 we must conclude that this simplest model does not satisfy the data. We 
look now for the least departure from the isotropic Gaussian model with non-zero axial dipole, 
that is consistent with the data. In order to do this it is worth restating the assumptions for this 
model and looking at whether it is possible to vary any of the basic parameters. The fundamen- 
tal assumption is that the spatial statistics of the present-day field are typical of the paleofield 
and may be used to determine the variance, a^, for the scaled Gauss coefficients of the non-dipole 
part of the field and also the non- axial parts of the dipole field. The magnitude of the axial part 
of the dipole field is also taken to be Gaussian, with the same variance, and a mean value given 
by the present day value for gi . There are three parameters that can be varied; 

1) a is computed from the present day field, and will depend on the level in the Earth at 
which the spectrum is assumed to be white, 

2) 7 ® the mean value for the axial dipole magnitude may not be the same as the present 
day value, and 

3) the variance for the dipole components of the field may be different from that for the 
non-dipole part. 

Since we have only directional data, varying a or yi will have the same effect, an increase in y ® 

corresponding to a decrease in a. Figure 10 shows the effect of varying yf on the model curves 
for B{i}, a(») and a(d) as a function of latitude. Here the parameters have been normalized so 

that yi° = 1.0 corresponds to the present day field value of 30p.T. In Figure 11 model curves are 
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plotted for various values of Vi, the standard deviation in each of the Gauss coefficients for the 
dipole part of the field (i.e., 7 ®, gi and hi). Oi is normalized in units of 7 ®, so that the giant 
Gaussian model whose pdfs are plotted in Figure 8 has oj =0.207 as determined from the spectrum 
of Figure 1 . It is readily seen that in order to make the model parameters a(») and o-(d) approxi- 
mate those for the data we must either increase yi (or equivalently decrease a) or decrease ctx. 
This also makes the zero crossing for B{t} closer to the right latitude, but makes the bias 
insufficiently negative at mid to high latitudes. There is little difference in the shapes of the 
curves obtained by varying these parameters, and the data is certainly not of good enough quality 
to enable us to pick between them. We must therefore look to other data to constrain these 
parameters. 

The variance of the scaled Gauss coefficients, and from these the variances in the field 
coefficients (see Section 1 and Appendix 2 ), are determined from the spectrum. Increasing the 
depth at which the spectrunii appears white, results in an increase in a and therefore a worse fit to 
the data. There is no evidence to suggest that the spectrum is white at a radius greater than the 
core-mantle boundary (in fact a number of authors have suggested it lies below there, e.g., Langel 
& Estes, 1982, find a level of 0.52a and a=36.7jjir, which results in a substantial increase in the 
variance in the field components). A decrease in the value of a is therefore not consistent with 
the known properties of the present field. 

Another source of information is the available paleomagnetic intensity data. McFadden & 
McElhinny (1982) have compiled the available data for the past 5Ma (with VGP latitudes of 
greater than 45°), and performed a statistical analysis of their virtual dipole moments (VDMs). 
The VDM is the equivalent dipole moment which would have produced the measured intensity at 
the calculated paleolatitude (assuming a dipolar field) of the sample. Working with VDMs 
enables these authors to combine data from sites at different latitudes. Unfortunately , there is 

no means of separating the gi component of the intensity from the non-axial dipole and non- 
dipole parts of the field. However, by taking random samples of the present day field, McElhinny 
& Senanayake (1982) concluded that the scatter in VDMs due to non-dipole components was 



- 18- 


Gaussian with standard deviations ranging between 12.5 and 20.7% of the mean. McFadden & 
McElhinny (1982) find a model for the VDMs consisting of nested distributions due to paleointen- 
sity errors, variation in the non-dipole parts of the field and a true dipole moment distribution. 
Their estimate for the preferred value of the true dipole moment corresponds to the peak (or 
mode) of the truncated Gaussian distribution used to model the true dipole moment, and is 
8.87 ±0.65 X 10** Am* (to within the 95% confidence limit). This corresponds quite closely to the 
value of 9.07 X 10** Am*, that they obtained from the arithmetic mean of the observed VDMs. 

The dipole moment p is related to the Gauss coefficients of degree 1 by 

P-— f(»f)“+(«.‘)‘+(k.M*r 

For the present day field p«7.91 x 10** Am*. If the and k}; terms are indeed Gaussian with 
zero mean, then their combined contribution to the true dipole moment will be that due to their 
scatter about the mean. Assuming, as before, that the present day field is typical, we can com- 
pare today’s dipole moment with McFadden & McElhinny’s estimate of the true dipole moment 

to see whether gi for the present day differs drastically from yi for the past 5Ma. The difference 

of about 10% suggests that this is not the case. We might be justified in increasing 7 ® by 10%, 
but hardly any more. Since the discussion above suggests that a is if anything underestimated 

(which would counteract the increase in 7 °) we choose to leave it unaltered. Figure 10 clearly 

shows that a 10% change in the value of 7 ® will not significantly improve the fit to the data. An 
increase of about 50% would be required in order to reduce the standard deviations in t and d to 
about the right level. This is completely inconsistent with what is known from both the paleoin- 
tensity data and the constraints on the present day geomagnetic spectrum. 

The constraints on 0*1 are much less rigorous. These come entirely from the present day 
field values. Our Gaussian model for the field assumes that the scaled dipole coefficients 7®, gi 
and hi , are drawn from the same distribution as the rest of the coefficients (except that 71^ has a 

non-zero mean). Support for this is based on the fact that and hi do not alter the 
Kolmogorov-Smirnov test for normality, when they are included in the distribution. However, 
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this may simply be chance; there are only two data involved, and there is considerable evidence 
for believing that the dipole part of the field is special. Why else would it dominate the non- 
dipole part so much of the time? Figure 11 shows that a value for o-i of 0.1 would provide a 
much more acceptable fit to the standard deviation data, than ux =0.207 obtained from the spec- 
trum. The present day values of gi: =0.0650ji® =0.1878^1® are entirely consistent with 

this, as both samples from the distribution lie within 2 standard deviations of the zero mesm. 
Decreasing cx still further [e.g. to 0.05, as shown in Figure 11) would enable us to constrain all 
the standard deviation data for the model to lie below those of the paleomagnetic data. However, 
this stretches the limits of what is plausible for the present day field, as it would require the 
present day sample of hi to lie at 3.764ax in the distribution. In a random sample from a Gaus- 
sian distribution this has a probability of less than one part in ten thousand of occurring, hardly 
what one would expect from a typical sample of the field. 

Reduction of the variance in the directional data, has been obtained at the expense of 
almost eliminating the bias in t. The zero crossing now occurs at about the right latitude, but the 
large negative values at mid-latitudes cannot be obtained from a dipole field with such low vari- 
ance. A different means of generating distributions with a significant bias in t is to include a 
zonal quadrupole term with a non-zero mean magnitude. 

The idea that the time averaged geomagnetic dipole does not correspond to that of a geocen- 
tric axial dipole (GAD) is not a new one: many paleomagnetists have proposed it (see e.g., Wil- 
son, 1971; Greer, Georgi & Lowrie, 1973; Wells, 1973; Georgi, 1974; Wilson & McElhinny, 1974; 
Merrill & McElhinny, 1977; Coupland & Van der Voo, 1980; Livermore, Vine & Smith, 1983, 
1984). Many of these authors have performed least squares fits for the spherical harmonic 
coefficients using paleopole data drawn from sites around the world. The resulting models may 
be rather strongly influenced by data distribution, whether non-zonal terms are included, and the 
degree to which the spherical harmonic expansion is carried. However, there appears to be a gen- 
eral consensus that for the past 5Ma the inclusion of a yj term offers a significant improvement 
over the GAD in fitting the available data. Livermore, Vine 8c Smith (1983) give a value of 
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7a = •057*1 agreement with that obtained by Merrill & McElhinny (1977). They also suggest 
that a value of 7s = .027* might be added for the last 5Ma, (although this could be due in p»t to 
data errors), and that hi may be significant in the average field for 0-5Ma. These results are 
substantially in agreement with a study of inclination anomalies by Lee & McElhinny (unpub- 
lished, but cited in Merrill & McElhinny (1983)), except that they suggest that the gS has a rev- 
ersing and non-reversing part. Lee & McElhinny also looked at the average values of gi and hi 
and found no evidence to indicate that they are significantly different from zero. They did not 
look at second order non-zonal terms in their analysis. 

Figure 12 shows the effect on the now familiar model parameters of varying 73 , the mean of 
the magnitude of 93 . All of these models have 7°= 1.0 and ai=0.1, as the evidence cited above 
suggests that these are reasonable values for these parameters. It is assumed also that gi is of 

the same sign as and reverses sign with gi . This is consistent with the paleomagnetic evidence 
that the inclination bias is negative for both normal and reverse fields (see e.g., Merrill & 

McElhinny, 1983, p.l85 figure 6.7). The best fit to the data of Lee (1983) is obtained with 73 =.06 
and is shown in Figure 13. Since our data set is somewhat sparse for our goal of determining the 
actual pdfs for directions at any given latitude, we mapped all the southern hemisphere data onto 

the upper hemisphere. Our modeling will therefore not indicate the necessity for a significant 73 
component, although it would be a simple matter to include such a term if it were necessary, and 
we had the data resolution to detect it. Similarly, we have made no attempt to include non-zonal 
components with non-zero mean in our model, because it is expected that by averaging data 
around latitude bands and over the 5Ma time interval non-zonal effects will average out. 

Figures 14 and 15 compare the cdfs for the preferred model (with 7®=1.0, o-i=0.1 and 

7s = .06) with the data distribution functions of Figure 6. Here, the solid lines represent the aver- 
age data cdfs over the latitude bands (0—15°, 15—30°, 30—45° and 45—70°) and the finely dashed 
lines give the model cdfs for the boundaries of the latitude strips. These should represent bounds 
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between which the data lie, provided the data are sufficiently good and the model is an adequate 
description of the statistics. For the modified inclination (Figure 14) the average model expected 
from combining data at the given sites is computed and is given by the heavy dashed line. It 
agrees well with the empirical cdfs, except in the 0—15^ strip for which either a larger quadrupole 
term or a lower dipole standard deviation would improve the fit somewhat. The declination data 
(Figure 15), as might be expected does somewhat less well. This is mostly because in both the 
15—30® and 45—70® strips the declinations are biased towards positive values. This could be due 
to bias arising from insufficient averaging in longitude in the sites used. Right handedness of 
VGP positions was also noted by Wilson (1972) and others since then. Lee (1983) finds no evi- 
dence for believing the time averaged g} and terms significantly different from zero, when all 
the sources of error in the data are taken into account. 

3. Conclusions 

In this study an attempt has been made to find a simple statistical model for the Gauss 
coefficients to describe the secular variation of the geomagnetic field. Working with the Gauss 
coefficients, instead of the directional data enables us to find a model that can simultaneously 
satisfy all the statistical requirements imposed by the vector field variation. Using the spatial 
spectrum of the present day field it is possible to derive a statistical model for the non- dipole part 
of the field. A comparison of the computed probability distribution functions for the modified 
inclination and declination data with paleomagnetic data spanning the past 5Ma indicates a 
number of features about the geomagnetic field statistics. (1) With the exception of the axial 
dipole and axial quadrupole terms each Gauss coefficient may be regarded as independently nor- 
mally distributed with zero mean and a standard deviation for the non-dipole terms consistent 
with a white source at the core surface. (2) The dipole part of the field cannot have the same sta- 
tistical distribution as the non-dipole part. The axial dipole has predominated and both the non- 
axial parts of the dipole and the magnitude of the axial part of the dipole have lower variance (by 
about a factor of 4) than if they belonged to the same giant Gaussian distribution as the non- 
dipole part of the field. (3) The hypothesis based on the present day field that all the non-dipole 
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Gauss coe£Bcients have zero mean is inconsistent with the paleodata. Although the inclination 
measurements are biased towards low values {relative to the GAD hypothesis) by the statistical 
variation in the Gauss coefficients, this bias alone is not sufficient to account for the inclination 
anomalies found in the paleodata; a quadrupole term with non-zero mean magnitude is required. 

A value of 73 =*067 ® for the past 5Ma provided a reasonably good fit to the paleodata; this agrees 
with typical values found by other workers (Livermore, Vine & Smith, 1983; Lee, 1983). 

The statistical distributions for the first 3 degree spherical harmonic coefficients are plotted 
in Figure 16 for this model. Note the quite narrow distribution of gi about its peak normal and 
reversed values, and how the variance drops off rapidly with increasing degree. The gS distribu- 
tion results from the sum of two normal distributions, centered at ±.06^i®. There is an implicit 

covariance between gi and gi , in that they are required to always have the same sign. Table 1 
lists some parameters associated with this model. This simple model for the Gauss coefficients 
also enables us to compute the expected probability density functions for the conventional field 
directions D and I at any latitude (Figure 17). The model is a complete description of the 
geomagnetic field in the sense that it is capable of predicting the expected statistics of any 
measurable feature of the field. For example, we could straightforwardly compute distribution 
functions for intensity measurements anywhere on the Earth; there is no decoupling of intensity 
and directional data in the spherical harmonic representation, making the model eminently 
testable by any or all of the available data. 

Some remaining questions clearly require further work. This study is based on a very limited 
data set of combined sedimentary and lava flow measurements. Information from lava flows is 
usually more suitable for statistical studies of secular variation, because there is less smoothing 
inherent in the recording process. A better data set requires finding all the available lava flow 
data for which individual flow mean values have been published and using these to obtain the 
empirical cdfs for a data set. This would require an extensive literature search since (to our 
knowledge) there is no existing compilation of such data (the published pole lists generally con- 
tain means from several sites rather than flow means for a given study). It remains unclear what 
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criteria should be used in selecting these data; internal consistency as represented by ags is clearly 
required, but it seems desirable not to exclude transitional data corresponding to low VGP lati- 
tudes, since these transitional fields are after all part of the secular variation. However, if all 
these data are included, we then encounter problems such as described in this paper with the Ice- 
landic data which appear to have too many transitional flows. This probably reflects the combi- 
nation of sampling two erratic processes, namely the reversal of the magnetic field and the vol- 
canic eruptions that generate the flows. A larger data set can only improve the situation. 

The acquisition of a large global data set could provide us with more accurate cdfs for i and 
d. This study has shown that it is possible to satisfy the data with a model of the type discussed. 
A better data set could give better estimates for the parameters of the model, perhaps justifying 
setting up a formal inversion procedure. 
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Appendix 1: Spherical Harmonic Representation of the Geomagnetic Field 

The method of analysis for global geomagnetic data was devised in 1838 by Gauss and 
involves representing the geomagnetic field of internal origin as the gradient of a potential T 
whose spatial variation is specified by an infinite sum of spherical harmonics. 

/ = lm=0 *’ 


B = -V'V 

i.e., 

o = n R = -iSl 

* r 90 rsin0 9<J> ' dr 

where <fi"^ and h/" are known as the Schmidt partially normalised Gauss coefficients, r, 0 and 4> 
are the usual spherical coordinates, and are the partially normalised Schmidt functions 
related to the associated Legendre polynomials, (Abramowitz & Stegun, 1965) by 


P,’" = P,„ for m=0 


pr = 




(/+m)! 


Pi^ for m>0 

The term "partial normalization" arises from the fact that within each degree / the average value 
of the square of P/”*over the surface of the sphere has the same value for all values of m, i.e. 


2ir 1 f j.'l r 

„ , |cosmd> , cosm d) 

/d<),/d(cos0)Pnco80)r ^ F,? (cose) 

0-1 v^y 


4'n 
2 / + 1 




The Gauss coefficients are commonly used to describe the field, since they are sufficient to deter- 
mine the field model at any point on the Earth’s surface. Successively higher degree and order 
terms correspond to components of the field with greater spatial variability; eg^ is the 
coefficient of the axial dipole term, g^ the axial quadrupole, etc. The zonal terms (m = 0) exhibit 
latitudinal spatial variation, and the non-zonal terms longitudinal variation. A more convenient 
mathematical representation of the field is in terms of the fully normalised complex spherical har- 
monics, y^”*(0,<})), with complex coefficients 
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^ = «i i (-rnrrne,<i>) 

/-Om--/ *■ 

The fully normalised spherical harmonics are related to the associated Legendre polynomials 
Pim(cose) by 

and 


2ir 1 

/J4./j(co89)y,i>'‘(e,*)r,"(e,<l>) = «r8,„' 

0 -1 

Although the Schmidt normalization is commonly used in geomagnetism, we will use this fully 
normalised representation where it simplifies the mathematics (see Appendix 2). The fully nor- 
malised tf”* may be related to the Schmidt coefficients by 



m>0 

m=0 


[(-ir{^T)‘ ^<0 

The Gauss coeffiicients are usually determined up to some finite degree by performing a least 
squares fit to the available field data (Langel, 1985, provides a review of methods and models). 
Downward continuation of these models enables us to make estimates of what the field looks like 
on the core-mantle boundary. Spherical harmonic fits can be highly susceptible to poor data dis- 
tribution, and this problem is exacerbated by downward continuation of the models. Shure, 
Parker and Backus (1982) developed a technique for finding the smoothest model, in a certain 
specified sense, that is consistent with the data. These models, known as harmonic spline models, 
to a large degree compensate for poor data distribution, by suppressing spatial variation in the 
field that is not required by the data. 


Appendix 2: Computation of Variances and Covariances for By, and 

As indicated in the previous section let us suppose that and A,”* are independent, normal 
random variables of zero mean for /a 2. Then y,'" and A" satisfy the following 
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E[gr\ = E\hr\ = 0 


E[grgr'\ = ^[w'l = 

1 = 0 for all 


w)iere <x^ is the variance of the normal distribution associated with degree /• First let us derive 
the variance of the non dipole part of the field components This is facilitated by 

using the fully normalised spherical harmonic representation for the potential. Making use of the 
relationships between the and the Schmidt coefficients, we find 


Tl 


urf+iw 




2 / + 1 
4tt 
2/ + 1 


O’/ 


Thus is the same regardless of whether or not m=0. Because |=0 for 1=1' 

or m = m ' the only non-zero terms possibly remaining are of the form (m 9^0), 

when gl^ and appear in both complex numbers. 

- 0 

We make use of the above expressions in deriving the field component variances and covariances. 
In addition it is useful to remember the Spherical Harmonic Addition Theorem (Jackson 1963, 
p.67): 

i rneo.Wi'ne.'H’ 

Letting f = fo yields 


i lyne,<l>)l' = ^ 

m--/ 4 tT 
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d 

Differentiating this equation yields some further sums which are useful. — yields 

30 


I n y’m 

2 = 0 
m« — / 


30 


Some further manipulation and differentiation of the addition theorem yields 

' 1 0 yr |2 _ /(f+i)(2/+i) 
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We start with the radial field. From our non-dipole field model £[Br]=0. Thus 
var(5,l = 

= 22S2(/+ 1)(' ' + 1) yi”X Yr' ' *] 

l r mm' 
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2 / + 1 
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,= 2 ^ (/ + 1)(2/ + 1)2 4ir 

= 

,t22l + Va^ 

a* c a 6a 

= .0753a* = (7.60)* (p,!)* 


Similarly £![fl(j]=0 and 


var 


f , -r-> dira,* _ , 3 y)”* ,, 

l^e] - ^2/ + i 5 30 

= 4ira* /(/ + 1)(2/ + 1) 

(2/ +!)*(/ + 1) 8ir 

2 « I 

2 i%(2l + iy a’ 

= .0262a* = (4.48)*(^t^)* 


Symmetry arguments clearly dictate that E[B^] = 0 and var(B^] = var(5g]. 
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Since and are derived from sums of Gaussian random variables, they are clearly 

Gaussian also. Their covariances may be readily computed from 


covlB,fl,| - = 2222(<+l)l'n-T^)'«IW"l 

I V mm' 

err 




1 = 2 

= 0 


Clearly by symmetry and similarly for E[B^B^], Thus the non-dipole field com- 

ponents B^j J00 and B^ are independent Gaussian variables with zero mean and standard devise 
tions which are equal for the two horizontal components and larger for the radial component. 


Appendix 3: Computation of the Distribution Functions for Declination and Inclina- 
tion 

Let us suppose that the pdf for the axial part of the dipole field, may be written as a sum 
of Gaussian kernel functions, Gy(rr), centered at positions Xj 

/x)(*) = <?,(*) 

i 

where 


G,(*) 


exp— % 



and 


This is a conventional technique of representing empirically an unknown pdf (see Silverman, 
1986); //)(x) is known as a kernel estimator. We will assume that Xy, the locations of the ker- 
nels, and Ay, their widths, are chosen a priori and hj = h — constant. Then the axial dipole 
contribution to the field components 5^ and 5,^ at any given site is distributed as 


F 

^ c,- f(a:.'-X,') ) 

/»(*') - 


2 



-29- 


where 


= k^z and = cos\ for 

Xj* == k^X^ and k^ = sin\ for 5, 

\ is the latitude at the site in question. 

In its simplest form (as is used in this work) //>(a:) is just the sum of two Gaussian distribu- 
tions, one centered on the value of the present day normal field, and the other on the reverse field 
(see the distribution function for Qi in Figure 16). Then if, as discussed in Section 1, the non- 
dipole part of the field has a Gaussian distribution, fsD{y)j independent of the dipole part 

of the field, the cumulative distribution function for the field components and B^ may be com- 
puted from 
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where S is the region satisfying x ' -I- y ^ and / includes the g / and h y random variation. 

Then / (p) the probability density function for the field component will be 




and will again be a sum of Gaussians 
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Note that this is easily modified to include a non-zero mean for the (and or any other higher 
order) term in the non-dipole distribution function. 

Now to compute the distribution function for tan/ we need first to find that for 
H = Letting 


Bq = /cos9 


= /sin0 


and 
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Yj = fcgXy = XjCOS\ 


we have 


/ 2TT 
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or for the pdf for H 
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The pdf for t = tan/ is given by 


/We (/cose-r,)^ 




/ta„/(0 = JWVr(«)//f(0 

0 

where is the pdf for the radial component of the field at i? = /ftan/, t.c. 
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and Zy = A^Xy = 2Xysin\ 


Some algebraic manipulation yields 
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/tan/ = exp(-|-)erfc(-^) 
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C = 


The 6 integral can be computed rapidly using the trapezoidal rule. The inclination pdf is then 

//(»■) = sec ^‘ / ta „/(0 


Similarly the declination distribution may be computed from 
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erfc(a;) is the complete error function (Abramowitz & Stegun, 1965, p.299) 

The pdfs for I and D are readily converted to those for the modified field directions i and d. 
Comparison with the data requires cdfs for d and i. These are obtained by using a cubic spline 
quadrature scheme to integrate the numerical values of the pdfs obtained on a uniform grid for d 
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and f\ Obtaining the cdfs also provides some reassurance about the accuracy of the above 
results, since we can check that they integrate to unity. 
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Table 1: Parameters describing preferred secular variation model. 
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Mean values for the spherical harmonic coefficients 
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30.0 1.8 0.0 0.0 0.0 


Standard deviation in surface field components (ti,T) 
5.39 5.39 9.68 


Standard deviation in non-dipole part of field components (|xT) 
4.48 4.48 7.6 
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Figure 2: Empirical cdf for the GSFC980 Gauss coefficients scaled 
in the manner described in the text. The dashed line shows the 
theoretical curve expected if the scaled coefficients correspond to a 
Gaussian distribution. 
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Figure 3: Empirical cumulative distribution functions for modified 
declination (d) and inclination (») data from Hawaii and Iceland. 
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Figure 4: Empirical cdfs for modified inclination, t, using the 
available individual site data from Lee (1983). All the data within 
each indicated latitude band are combined to generate the cdf. 
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Figure 5: Bias in the modified inclination for the data grouped as in 
Figure 4 (lower part of figure). Data from the two smallest data 
groups are marked by a star to indicate their lesser reliability. Upper 
part shows the standard deviation in d (triangles) and i (squares) for 
the same data. Error bars are one standard error computed using the 
bootstrap technique described in the text. 
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Figure 6: Empirical cdfs for t from individual site data of Lee 
(1983), with data grouped into wider latitude bands than in Figure 4. 
Note the smoother cdfs obtained. 
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Figure 7: Bias for the modified inclination and standard deviations 
for d and i for the data as grouped in Figure 6. 
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Figure 8: Probability density functions for d and i at a variety of 
latitudes, assuming the giant Gaussian model for the secular 
variation, with a mean axial dipole corresponding to the present day 
value. 
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Figure 9; Bias and standard deviations as a function of latitude for 
I and d for the model (solid line) whose pdfs are shown in Figure 8 
and the data of Figure 4 (symbols). 
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Figure 10; Effect of varying on the bias and standard deviation 
for the resulting i and d distributions. 
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Figure 11s Effect of varying <Ti, the dipole standard deviation, on 
the bias and standard deviation for t and d. 








Figure 12: Bias and standard deviation curves resulting for various 
values of 
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Figure 13; Bias and standard deviations as a function of latitude for 
the preferred model described in the text. Parameters for the data 
compiled by Lee (1983) are again shown as open symbols. 
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Figure 14: Cdfs for i data of Lee (1983) grouped as in Figure 6 
(solid line) compared with the predictions of the preferred model. 
Long dashed lines indicate the expected average cdf from combining 
data at the given sites. The finely dashed lines indicate the region 
within which the observed cdf should lie. 
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Figure 15: Cdfs for d data of Lee (1983) grouped as in Figure 6 
(solid line) compared with the predictions of the preferred model. 
The finely dashed lines indicate the region within which the observed 
cdf should lie. 
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Figure 16: Statistical distribution for the first 3 degree spherical 
harmonic coefRcients for the model. Parameters for these 
distributions are given in Table 1. 
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Figure 17: Probability distribution functions for declination, D, and 
inclination, / at a variety of latitudes for the preferred model. 




